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Abstract 

Hopf bifurcations in fast-slow systems of ordinary differential equations can be 
associated with surprising rapid growth of periodic orbits. This process is referred 
to as canard explosion. The key step in locating a canard explosion is to calculate 
the location of a special trajectory, called a maximal canard, in parameter space. A 
first-order asymptotic expans ion of this location was found by Krupa and Szmolyan 



Krupa and Szmolyanl . l2001bl jal)d] in the framework of a "canard point" -normal- form 



for systems with one fast and one slow variable. We show how to compute the coef- 
ficient in this expansion using the first Lyapunov coefficient at the Hopf bifurcation 
thereby avoiding use of this normal form. Our results connect the theory of canard ex- 
plosions with existing numerical software, enabling easier calculations of where canard 
explosions occur. 

1 Introduction 



Our framework in this paper is the theory of fast-slow ordinary differential equations (ODEs): 



ex 



CLT 



dy 



y = -j^ = 9ix,y,X,e) 



where (x, y) G x R", A G M is viewed as a parameter and e is sufficiently small, i.e. 
< e < 1. The functions / : M"* x M" x M x R ^ M"* and ^ : x M" x M x M ^ R'^ 
are assumed to be at least in this paper. The variables x are fast and the variables y are 
slow. An in troduction to the theory of fast-slow systems fro m the geometric viewpoint can 
be foun d in Arnold . 1994 : Jones . 1995 : Guckenheimerl. 2002 ] . asymptotic methods are devel- 
oped in Mishchenko and Rozovl. 1980| : Grasman . 1987 ] and ideas from nonstandard analysis 



are considered in 
methods here. 



Diener and Dienen . Il995| . We will only use geometric and asymptotic 



In the singular limit e — )■ the system ([T]) becomes a differential-algebraic equation. The 
algebraic constraint defines the critical manifold: 

Co = {(x,y)GM'"xR":/(x,i/,A,0) = 0} 
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For a point p G Cq we say that Cq is normally hyperbolic at p if all the eigenvalues of the 
m X m matrix D^fip) have non-zero real parts. A normally hyperbolic subset of Cq is an 
actual manifold and we can locally parametrize it by a map ip{y) = x. This yields the slow 
subsystem (or reduced flow) ij = g{'ilj{y),y, A, 0) defined on Cq. 

Changing in ([1]) from the slow time scale r to the fast time scale t = r/e yields: 

x' = ^ = /(x,?/, A,e) (2) 

dy /XX 
y = = ^9{x,y,X,e) 

Taking the singular limit e — )■ in ([2]) gives the fast subsystem (or layer equations) x' = 
f{x, y, A, 0) with the slow variables y acting as parameters. A point p G Cq is an equilibrium 
of the fast subsystem. We call a subset S' C Cq an attracting critical manifold if all points p 
on it are stable equilibria of the fast subsystem i.e. all eigenvalues of D^fip) have negative 
real parts. The subset S* C Cq is called a repelling critical manifold if for all p G S* at least 
one eigenvalue of D^f{p) has positive real part. 



Fenichel's Theorem [Fenichell . Il979| states that normally hyperbolic critical manifolds 
perturb to invariant slow manifolds C^. A slow manifold is 0(e) distance away from Cq. 
The flow on the (locally) invariant manifold converges to the slow subsystem on the crit- 
ical manifold as e — )■ 0. Slow manifolds are usually not unique for a fixed value of e = eo but 
lie at a distance 0{e~''^'^") away from each other for some k > 0; nevertheless we shall refer 
to "the slow manifold" associated to subset of the a critical manifold with the possibility of 
an exponentially small error being understood. 

Suppose the critical manifold can be divided into two subsets Sa and Sj. where Sa is 
attracting and Sr is repelling so that Cq = SaU LU Sr- Here L denotes the part of Cq that 
is not normally hyperbolic. We assume that for p & L the matrix Dxf{p) has a single zero 
eigenvalue with right and left eigenvectors v and w and that w ■ Dxx{p){v^ v) and w ■ Dyf{p) 
are non-zero. In this case points in L are called fold points. We can use the flow of ([2]) to 
extend the associated slow manifolds Sa,e and Sr,e but the extensions might not be normally 
hyperbolic. The key definition used in this paper is that a trajectory 7 in the intersection 
of Sa,e and Sr^e IS Called a maximal canard; note that this definition requires the extensions 
of the slow manifolds under the flow. Observe that 7 C Sr^e despite the fact that Sr^e is 
repelling in the fast directions. 

We are interested in the case when a fast-slow system undergoes a Hopf bifurcation and 
a maximal canard is formed close to this bifurcation. The periodic orbits resulting from the 
Hopf bifurcation grow rapidly in a A- interval of width 0{e~^/'^) for some > 0. The rapid 
orbit growth is usually referred to as canard explosion and the bifurcation scenario is called 
singular Hopf bifurcation. 

The paper is organized as follows. In Section [2] we describe re sults on singular Hopf 



bifurcation and canard explosion obtained by Krupa and Szmolyan [Krupa and Szmolyanl . 
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2001bl | . In Section |3] we clarify the different definitions of the first Lyapunov coefficient of 
a Hopf bifurcation. In Section H] we present the main results on the relation between the 
location of the maximal canard and the first Lyapunov coefficient. We describe which terms 
will contribute to the first order approximation using a rescaled Hopf bifurcation normal 
form. Then we show explicitly how to compute a first order approximation to the location of 
the maximal canard avoiding additional center manifold reduction and normal form trans- 
formations. In Section we locate the maximal canards in two examples: a two-dimensional 
version of van der Pol's equation and a three-dimensional version of the FitzHugh-Nagumo 
equation. 



Note that we do not give a detailed description of dynamics associated to a singular Hopf 
bifur c ation and refer the reader to the previous ex t ensive literature e.g. Baer and Erneuxl . 



19861 . 11992I : iBraaksmal . Il998l : iKrupa and Szmolyad . l2001bl : iGuckenheimerl . 12008 



2 Canard Explosion 



We describe the main r esults about canard explosion in fast-slow systems with one fast and 
one slow variable from Krupa and Szmolyaru . l2001bj . Consider a planar fast-slow system of 
the form 



X 

y' 



A,e) 
(^9{x,y,X,e) 



(3) 



where f,geC'^ 



'-, M) for k > 3, AGMisa parameter and < e ^ 1. Denote the critical 
manifold of ([3]) by Cq. We assume that Co is locally parabolic with a minimum at the origin 
{x,y) = (0,0) independent of A so that (0,0) is a fold point; more precisely 



/(0, 0, A, 0) = 0, UO, 0, A, 0) = 0, A.,.(0, 0, A, 0) ^ 0, /,(0, 0, A, 0) ^ 



(4) 



In addition, we assume that g{0, 0, A 7^ 0, 0) 7^ 0; under these conditions the fold point at 
the origin is generic for A 7^ 0. We assume without loss of generality that /xx(0, 0, A, 0) > 
so that Co is locally a parabola with a minimum at the origin. Using and the implicit 
function theorem we have that Co is the graph of a function y = (j){x) for :[/—)■ R where 
f/ is a sufficiently small neighbourhood of x = 0. Assume that Co splits into an attracting 
and a repelling curve C = Q U {(0, 0)} U where 



Q = {x<0,f,<0}nCo, Cr 

The situation is shown in Figure [U^a). Differentiating y 
we get that the slow flow on Co is defined by 

dx . (7(x, 0(x), A, 0) 
dr 0'(x) 



{x>o,/,>o}nCo 

(/)(x) with respect to r 



te 



Note that the slow flow is singular for A 7^ at (0,0) since 0'(O) = and (7(0, 0, A, 0) 7^ 0. 
Assume that at A = we have a non-degenerate canard point (see Figure [It^b)) so that in 
addition to the fold conditions we have 



(7(0, 0, 0, 0) = 0, (7,(0, 0, 0, 0) ^ 0, (7a(0, 0, 0, 0) 7^ 
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Figure 1: (a) A generic fold for A 7^ 0. (b) A nondegenerate canard point for A = 0. The 
slow flow is indicated by single and the fast flow by double arrows. 

Therefore the slow flow is well-defined at (0, 0) for A = and we assume without loss of 
generality that x > in this case. Near a rion-dege nerate canard point ([3]) can be transformed 
into a normal form Krupa and Szmolyanl . 2001c| : 



x' = -yhi{x,y,X,e) + x'^h^x,y,\,e) + ehs{x,y,\,e) 

y' = (^ixh4ix,y,X,e) - Xh{x,y,X,e) +yhQ{x,y,X,e)) (5) 

where the functions hi are given by: 

hix,y,X,e) = 0{x,y,X,e) 

hj{x,y,X,e) = 1 + 0(x, A, e), j = 1,2,4,5 

We define several computable constants, abbreviating (0, 0, 0, 0) = in the definitions: 

ai = {h)x{0), a2 = {hi)^{0), as = (/i2)z(0), 04 = {h^)x{0), 05 = (/i6)x'(0) 

Note that all a, for i = 1,2, 3, 4, 5 only depend on partial derivatives with respect to x. Next 
we define another constant: 

A = — a2 + 803 — 2a4 — 2a5 



Theorem 2.1. ( jKrupa and Szmolyan . \2001a] ) For < e < eo, |A| < Aq and eo > 0, Aq > 



sufficiently small and under the previous assumptions in this section there exists a unique 
equilibrium point p for in a neighbourhood of {x, y) = (0, 0) . The equilibrium p undergoes 
a Hopf bifurcation at Xh with 

XH = -^^e + Oie'/') (6) 
The slow manifolds C^^i and C^^r intersect/coincide in a maximal canard at Ac for 

^c = -(^^ + i)e + 0{e'/') (7) 



The equilibrium p is stable for X < Xh and unstable for X > Xh- The Hopf bifurcation is 
non- degenerate for A 0, supercritical for ^ < and subcritical for A > 0. 
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Remark: The asymptotic expansions for Xh and Ac are asymptotic series with asymptotic 
sequence {e'^/^j^Q and Theorem 12.11 imphes that the first two coefficients of the expansion 
are zero and the third coefficient can be computed exphcitly. 



Note that currer itly no standard bifurcation softw are such as AUTO [Doedeh Champneys. Dercole. Fairs 



20071 or MatCont [Govaerts and Kuznetsovi . |2008| | computes the constants and A auto- 



matically. Nevertheless bifurcation software can detect Hopf bifurcations so that given a 
fixed e we can approximate Xh numerically. Hence the numerical problem that remains is 
to compute A since 

A 



Xh — Xr 



To simplify the notation we define K = A/8. If we know K we can easily approximate the 
location of the maxir nal canard by An = — Ke + 0(e^^^). The maximal canard organizes 
the canard explosion Krupa and Szmolyanl . l2001bl | and indicates where the rapid amplitude 
growth of the small orbits generated in the Hopf bifurcation occurs. Our goal is to avoid 
any additional normal form transformations and center manifold reductions to compute K. 
The key point to achieve this is to observe that K is just a rescaled version of "the" first 
Lyapunov coefficient of the Hopf bifurcation at Xh- 



3 The First Lyapunov Coefficient 

We review and clarify the interpretation, computation and conventions associated with the 
first Lyapunov coefficient of a Hopf bifurcation. Consider a general N-dimensional ODE at a 
non-degenerate Hopf bifurcation point. We assume that the equilibrium has been translated 
to the origin so that 

z' = Mz + F{z), for zeR^ (8) 
with F{z) = 0(||2;|p) and M = {rriij). Taylor expanding F yields 

z' = Mz + Ib{z, z) + ]:C{z, z, z) 
where the multilinear functions B and C are given by: 



N 



Bi{u,v) 
Ci{u,v,w) 



h ^^^^^^ 



UjVk 



N 

j,k,i=i 



UjVkWi 



The matrix M has eigenvalues Ai,2 = ±«wo for uq > 0. Let q G C be the eigenvector of Ai 
and p G the corresponding eigenvector of the transpose M'^ i.e. 



Mq = iuoq, Mq = —iuoq, M^p = —iuop, Mp = iuop 

where the the overbar denotes componentwise complex conjugation. We can always normal- 
ize p so that the standard complex inner product with q satisfies p^q = J2f=iPj1j — 1- The 
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first Lyapunov coefficient of the Hopf bifurcation can then be defined by ( jKuznetsov . 2004 
p.180): 



Ku 



2Ur 



{p^C{q, q, q) - 2p^ B{q, L-'B{q, q)) + p^B{q, {2iuoIn - M)-'B{q, q))) (9) 



In the case of a two- dimensional vecto r field F = {F^, F"^) the formula (Q can be expressed 
in the simpler form ( [Kuznetsov . 2004 1. p. 98): 



Ku 



2ui 



(10) 



where 



g2o = p^B{q,q), gu=P^B{q,q), g2i = P^C{q,q,q) 

It is important to note that /f"" is not uniquely defined until we choose a normalization of the 
eigenvector q. We adopt the convention using unit norm q^q = 1. A slight modification of the 
formula is used to evaluate the Lyapunov coeffi cient /{^*" numerically in the bifurcation 
software MatCont Govaerts and Kuznetsov . 2008 1. Using the current MatCont conventioE0 
we note that 

Other expressions for the first Lyapunov coefficient can be found in the literature. We 
consider only the planar case using simpler notation {zi,Z2) = {x,y): 

x' \ F\x,y) 
y' \F\x,y) 



--: M 



+ 



fix,y) 
9{x,y) 



Then another convention for li is ( |Chow. Li. and Wang . 1994| . p. 211): 



CLW 



^12 . 2 



^alifxxx + Qxxy) + 2m22{fxxy + dxyy) ~ ^2l{fxyy + Qyyy)] 



—mi2m22{fxx ~ fxxdxy ~ fxyQxx ~ QxxQyy ~ "^Qxy) 
-m2im22{gly - Qyyfxy " Qxyfyy " fxxfyy " '^fly) 
'^fTT'l2{fxx9xx + 9xx9xy) "^21 (/j/J/fi's/y ~^ fxyfyy) 



(12) 



(uJq + 3^7122) ( /ax /xy dxydyy)] 



where all evaluations in f lT2|) are at {x,y) = (0,0). Next, assume that we have apphed a 
preliminary linear coordinate change 



where A^ : 



to the system (ITTl) to transform M into Jordan normal form. Then we look at: 



u 
v' 
























+ 



g{u,v) 

f*{u,v) 
g*{u,v) 



(13) 



^MatCont version 2.5.1 - December 2008 
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In this case the Lyapunov coefficient formula simphfies ( jCuckenheimer and Holmes , 1983 
p.152): 



GH 



-j^g xxx ~^ f xyy ~^ 9xxy 



+ 9yyy\ + 



If* (f* f* 
U xy\J XX ' J y'l 



yy' 



9xyi.9xx ~^ 9yy) f xx9xx ~^ fyy9yy\ (-^^) 

Note that the linear transformation N is not unique. We adopt the convention that 



2Re(gi) 
2Re(g2) 



-2Im(gi) 
-2Im(g2) 



where q= (gi, ^2) is the normalized eigen vector of the linearization L that satisfies Lq = iuoq. 
Another common definition for (fT3ll is ( jPerko . 2001], p. 353): 



Pe 



37r 

—(If* f* + f* a* - f* 0* - a* a* - 0* 0* + f* f* ] 

^'-•> xy J yy ' •'yy^yy Jxxjxx Jxyjxx Jxyjyy ' JxyJxxl 



(15) 



~^^0 \9yyy ~^ f xxx ~^ f xyy ~^ 9xxy\ ) 

The Hopf bifurcation theorem holds for any version of li as only the sign is relevant in this 
case: 



Theorem 3.1. (see e.g. iGuckenheimer and Holrnea . \l98di : \Kuznetsou . \2004j ) A non-degenerate 
Hopf bifurcation of ([8]) is supercritical if li < and subcritical if li > 0. 



Since we need not only a qualitative result such as Theorem I3.H but a quantitative one 
relating the Lyapunov coefficient to canard explosion, it is necessary to distinguish between 
the different conventions we reviewed above. 



4 Relating li and K 



Krupa and Szmolyan consider a blow-up Krupa and Szmolyan . 2001b 313] of ([5]) given by 
$ : 5''^ X / — in a particular chart K2: 



X 



rx2, 



y = ry2, 



X = rX 



2, 



r €2 



(16) 



where r G / C M and (^2, 1/25 ^2, €2) G S^. Using flTB]) the resulting vector field can be 
desingularized by dividing it by y/e. We shall not discuss the details of the blow-up approach 
and just note that this transformation and the following desingularization are simply a 
rescaling of the vector field given by: 

e-'/'X, 



X2 



2/2 



Xo 



t2 = e'/H 



(17) 



Using the formula from Chow, Li and Wang [Chow. Li. and Wangl . Il994| in the rescaled 
version of ([5]) Krupa and Szmolyan get the following result: 

Proposition 4.1. In the coordinates ( ITTl) the first Lyapunov coefficient li has asympototic 
expansion: 

JCLW = K^e + 0{e) (18) 
where the overbar indicates the first Lyapunov coefficient in coordinates given by ( IT71) . 
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First, we want to explain in more detail which terms in the vector field ([5]) contribute to 
the leading order coefficient K. The main problem is that the Lyapunov coefficient is often 
calculated after an e-dependent rescaling, such as fll7l) . has been carried out. This can lead to 
rather unexpected eff ects in which terms c ontribute to the Lyapunov coefficient, as pointed 



out by Guckenheimer (Guckenheimerl . l2008l | in the context of singular Hopf bifurcation in 



To understand how the rescaling f ll7p affects the Lyapunov coefficient we consider the 
Hopf normal form case. We start with a planar vector field with linear part in Jordan form 
( !T3|l . Assume that the equilibrium is at the origin (x, ?/) = and Hopf bifurcation occurs for 
A = 0. Applying the rescaling (fT7|) we get: 

dx2/dt2\^f -uo/^/e \ f X2 \ f -J*{^/ex2,ey2) \ , . 

dy2/dt2 J V ^o/V~e j{y2 Jj,g*{V~ex2,ey2) J ^'^^ 

In a fast-slow system with singular Hoof bifurcation we know that g*{., .) = e(. . .) and that 
Uq = 0{y/e). Setting k^j = wo/v^ the Lyapunov coefficient can be computed to leading order 
by (HI: 

^T"" = j- (0, 0)K,.,(0, 0) + f:^jO, 0)] + /c./.,.,.,(0, 0)) + 0(e) (20) 

Equation ( 120|) explains the leading-order behaviour more clearly and shows that due to the 
rescaling certain derivative terms in the Lyapunov coefficient for a singular Hopf bifurcation 
are non-leading terms with respect to e — )■ 0. The point is that the rescaling modifies the 
order with respect to e of the linear and nonlinear terms. Also, applying the chain rule to 
the nonlinear terms to calculate the necessary derivatives can affect which terms contribute. 

To make Proposition (14. ip more useful in an applied framework we have computed all 
the different versions of the Lyapunov coefficient defined in Section ([3]) up to leading order 



for eq 



Inc 



natio n ([5]) in original non-rescaled coordinates. The computer algebra system Maple 



20081 was used in this case: 



^f^ = ^ + 0(7^) 

/f^ = K + Oie) (21) 
^CLW ^ K + Oie) 

= |f^ + 0(eM) 

64:UJ0 

Using the results ( 12T]1 we now have a direct strategy how to analyze a canard explosion 
generated in a singular Hopf bifurcation. 

1. Compute the location of the Hopf bifurcation. This gives Xh- 

2. Find the first Lyapunov coefficient at the Hopf bifurcation, e.g. we get Zf^"" ~ 4KuJo/^/e. 
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3. Compute the location of the maximal canard, and hence the canard explosion, by 
Ac ~ — Ke. For example, using MatCont we would get 

iMC 

Xc^Xh- -^e'/' (22) 

Observe that the previous calculation may require calculating the eigenvalues at the Hopf 
bifurcation to determine Uq but does not require any center manifold calculations nor addi- 
tional normal form transformations; these have basically been encoded in the calculation of 
the Lyapunov coefficient. 



5 Examples 



The fir st example is a version of van der Pol's equation |der Poll . ll920lll926l : lKrupa and Szmolyanl . 
2001b| given by: 



X 



y-x - 



X 



y' = e(A - x) 

We have to reverse time t — )■ — t to satisfy the assumptions of Section ([2]). This gives: 



X 



X 



X 

e(x — A) 



(23) 



The critical manifold is given by Cq = {y = x"^ + x^/3} with two fold points at (0,0) and 
(—2,4/3). The fold points split the critical manifold into three normally hyperbolic parts: 

Ci = Con{x< -2}, Cm = Con{-2<x < 0}, c, = Co n {0 < x} 

We only study the fold point at the origin which becomes a canard point for A = 0. The 
unique equilibrium point p = (xe(A), ye(A)) of (!23l) lies on Cq and satisfies Xe(A) = A. It is 
easy to check that subcritical Hopf bifurcation occurs for A = Xh = 0. Matching terms in 
and the normal form we find: 



hi = /i4 







Therefore K = 1/8 and we find analytically that the location of the maximal canard repre- 
senting the intersection of Cm,e and (7^,^ is 

Ae = -(l/8)e + 0(e3/2) (24) 
A numerical continuation calculation using a bifurcation software tool - we used MatCont 



Govaerts and Kuznetsovl . |2008| - gives that the first Lyapunov coefficient for the Hopf bi- 



furcation at A = Ajy = for e = 0.05 is 



MC 



0.4762 
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An easy calculatioij^ yields that ujq ~ 0.2236. Using ( 122|) we compare this to the result in 
equation f l24p with e = 0.05. Dropping higher-order terms we have: 

Ac(analytical) = —0.0063, Ac(numerical using li) = —0.0060 (25) 

The coincidence of the values of the location of the maximal canard is already quite good 
but this is expected since we have only compared the asymptotic formula to the Lyapunov 
coefficient formula derived from it which was evaluated numerically using continuation. A 
simple direct test to compare fl25|) to the location of the maximal canard is to use continua- 
tion of periodic orbits from the Hopf bifurcation point. The results are shown in Figure |2j 




Figure 2: Continuation of periodic orbits emanating from the Hopf bifurcation at A = 0. 
The parameter values for the red orbits are A = -0.001, -0.0025, -0.004, -0.005, -0.006, 
—0.0065 and for all the green orbits the parameter value is A ~ —0.006509 indicating a 
canard explosion near this parameter value. 



Remark: Depending on the bifurcation software used, direct continuation of periodic 
orbits can fail for small values of e. In this case spe cial methods are needed t o con- 
tinue periodic orbits having canard segments; see e.g. Guckenheimer and LaMarl . 120071 : 



Guckenheimer and Kuehij . l2009bl : iDesroches. Krauskopf. and Osingal . |2010|. Note that lo- 



cating Hopf bifurcations and calculating Lyapunov coefficients works well even for very small 
values of e as we require only local algebraic calculations. 



We conclude from Figure [2] that our estimates in (1251) are very good indicators to deter- 
mine where the canard explosion exists since they are already decent for a relatively large 
e = 0.05. In many standard fast-slow systems values of e < 0.01 are commonly considered. 



^Using MatCont 2.5.1. we can modify the file /matcont2.5. 1/MultilinearForms/nf _H.m to return the 
variable omega= ojq or to return = if^'-^ /loq. 
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In higher dimensions the analytical calculations will be very di fficult to carry out. As a 

second example consider a version of the FitzHugh-Nagumo equation Champneys. Kirk. Knobloch. Oldem^ 



20071 : iGuckenheimer and Kuehnl . l2009al . 12010 



Ob n 



y 



X2 
1 . 



SX2 — Xi{xi — 1)(0.1 — Xi) + y — I) 



(26) 



(xi - y) 



Equation (1261) has two fast and one slow variable and a 
p. For a detailed fast-slow system a r ialysi s describing the 

We only note 



where / and s are parameters, 
unique equilibrium point p{I) = 

bifurcations we refer the reader to Guckenheimer and Kuehnl . l2009al . 12010 
that the critical manifold is cubic curve given by 

Co = {{xi, X2, y) G M'^ : X2 = and y = Xi{xi — 1)(0.1 — Xi) + /} 

It is normally hyperbolic away from two fold points xi^± given by the local minimum and 
maximum of the cubic. The equilibrium p passes through the fold points under param- 
eter variation; 0(e) away from these p oints the equilibrium undergoes Hopf bifurcation 
Guckenheimer and Kuehnl . l2009al . |2010|. We shall just compute a particular case applying 
our result fl2T]) . We fix s = 1.37 and observe that / plays the same role as A in our previous 
calculations. We calculated the location of the maximal canard for several values of e. The 
results are shown in Table [H 



e 


Ic 


Ic{Lyapunov) 


10-2 


^ 0.0582046 


^ 0.06308 


5 • 10-3 


^ 0.0545535 


^ 0.05629 


10-3 


^ 0.0517585 


^ 0.05196 


5 ■ 10-4 


^ 0.0514108 


^ 0.05150 



Table 1: Comparison between the actual location of the maximal canard (canard explo- 
sion) Ic and the first-order approximation Ic{Lyapunov) computed using the ffist Lyapunov 
coefficient at the Hopf bifurcation. 



The second column of Table Hjshows the actual location Ic of the r aaximal canard (canard 

explo sion) obtained from continuation of periodic orbits using AUTO Doedel. Champneys. Dercole. Fairgri 
2007|. The third column shows the approximation obtained by using the first Lyapunov coef- 
ficien t. The Lyapunov coefficient has been computed using MatCont [Govaerts and Kuznetsovl . 



2008| . The error E{e) of this calculation is of order 0{e^^'^) as expected from ()18p . Obviously 



the approximation improves for smaller values of e. 



In the case of e = 0.01 it has been shown in [Guckenheimer and Kuehnl . I2009al . l2010l | that 
there is an intricate bifurcation scenario involving homoclinic orbits in a parameter interval 
near Ic- The first order approximation of the maximal canard is not sufficient to relate 
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it to the homoclinic bifurcation. This shows that the magnitude of e and other relevant 
bifurcations in the system have to be taken into account carefully when applying the results 
we presented here. 

6 Discussion 

We have investigated the relation between the first Lyapunov coefficient at a singular Hopf 
bifurcation and the associated maximal canard orbit. The major result is that no additional 
algorithms are needed to compute a first order approximation to the location of the maximal 
canard. Standard bifurcation software packages compute the Lyapunov coefficient and our 
results can be used to approximate the maximal canard location from this numerical calcu- 
lation. 

We also pointed out that there is no "standard definition" of the first Lyapunov coef- 
ficient of a Hopf bifurcation. This is not surprising since classical qualitative bifurcation 
theory only requires the sign of the Lyapunov coefficient. We hope that the comparison in 
Section [3] will help the reader to adapt their own numerical algorithms and software packages 
to support the calculation of maximal canard locations. 

Open questions which we leave for future work include the extensions to multiple slow 
variables, higher-order asymptotic expansions and the relation between the Lyapunov coef- 
ficient and blow-up transformations. 
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